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ABSTRACT 

Using iV-body simulations, we studied the detailed evolution of central stellar velocity dis- 
persion, (7*, during dissipationless binary mergers of galaxies. Stellar velocity dispersion was 
measured using the common mass- weighting method as well as a flux-weighting method designed 
to simulate the technique used by observers. A toy model for dust attenuation was introduced 
in order to study the effect of dust attenuation on measurements of cr* . We found that there are 
three principal stages in the evolution of a,, in such mergers: oscillation, phase mixing, and dy- 
namical equilibrium. During the oscillation stage, cr* undergoes damped oscillations of increasing 
frequency. The oscillation stage is followed by a phase mixing stage during which the amplitude 
of the variations in cr* is smaller and more chaotic than in the oscillation stage. Upon reaching 
dynamical equilibrium, cr* assumes a stable value. We used our data regarding the evolution 
of during mergers to characterize the scatter inherent in making measurements of cr* in non- 
quiescent systems. In particular, we found that o; does not fall below 70% nor exceed 200% of 
its final, quiescent value during a merger and that a random measurement of cr^ in such a system 
is much more likely to fall near the equilibrium value than near an extremum. Our toy model of 
dust attenuation suggested that dust can systematically reduce observational measurements of 
oi and increase the scatter in cr* measurements. 

Subject headings: Galaxies: evolution, Galaxies: interactions, Galaxies: kinematics and dynamics. Meth- 
ods: numerical 



1. INTRODUCTION 

The central stellar velocity dispersion, cr,, of a 
galaxy is a key observable quantity in galactic as- 
tronomy. Its importance is primarily due to the 
fact that (T» can be used as a proxy for the gravita- 
tional pot ential when in terpreted using the Jeans 
equations (|Jeans 1119151) or the Virial theorem. In 
particular, there are two important galaxy scal- 
ing relations involving cr* which are of great im- 
portance to the study of galaxy formation and 
evolution: th e fundamental plane (FP) of ellip 
tical galaxies (Dressier et al. Il987t iDavies et al. 



1987t Bender et al. 1992) and the Mbh-Q-* rela- 



tion ( Ferrarese fc Merritt I jgooot iGebhardt et al. 
2OOOI : iTremaine et al. 1120031 



The FP can be expressed in various forms, most 
commonly L oa aj So where L is the luminos- 



ity of the galaxy and Sc is the average surface 
brightness within the half-light radius (Re), or 
alternatively Rc oc cr^'^/"'^'^ where the intensity, 
/c = I{Rc)- The relation is significant because it 
can be used to estimate the distance to an elliptical 
galaxy, but perhaps more importantly, because its 
existence yields clues to how elliptical galaxies are 
formed. Some studies have suggested that ongo- 
ing mergers and recent merger remnants do not fit 
well onto the FP relation, while other studies show 
that late-stage mergers (i.e., mergers containing a 
single nucleus) can fit near or directly on the FP 
relation if several complicating factors are prop- 
erly taken in to account (fo r detailed di scussions. 



see Shi er fc Fischer 



Dasvra et al 



Rot hb erg et al. 



2006|: 
2OI1I 



1998t iBender et al. 



Rothberg fc Joseph I 12006 



and references therein). 
The two primary complicating factors that arise 
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when studying the FP relation in ongoing and 
recent mergers are (1) the presence of dust ob- 
scuration and (2) enhanced star formation, which 
temporarily increases the relative abundance of O 
and B stars. Therefore, to gain a better under- 
standing of the FP relation — especially in merger 
remnants — stellar dynamics must be disentangled 
from these complicating factors. 

The Mbh'C* relation is a tight relation between 
the mass of the central black hole (Mbh) and cr*, 
of the form Mbh oc a". There is evidence that 
black hole (BH) activity may play a role in galaxy 
formation and evolution, e.g., by regulating star 
formation through winds and outflows. To un- 
derstand galaxy formation, we therefore need to 
understand how super-massive BHs and galaxies 
co-evolve. The MBH~cr* relation is one of the 
most important tools available in this endeavor. 
In order to fully understand the origin and im- 
plications of the A/bh~ct* relation, we must deter- 
mine how the relation evolves with cosmological 
time (i.e., redshift). Active galactic nuclei (AGN) 
have become instrumental in studying the evolu- 
tion of the A/bh^ct* relation with redshift because 
AGN hosts are the only galaxies for which we can 
measure A/bh in the non-local universe. Unfor- 
tunately, limiting the study of A/bh-c* to AGN 
hosts may introduce a bias or significantly increase 
the observed scatter. This is because AGN activ- 
ity has often been linked with gala xy merger ac- 
tivity (jCanalizo fc Stockton I l20Ql[ ). Many mod- 
ern numerical simulations that successfully repro- 
duce present-day properties of early-type galaxies, 
such as the color-magnitude or the AfBH^cr* rela- 
tions, presuppose that AGN acti vity is triggered 



by mergers of gas-rich galaxi es (jHopkins et al 



200d [Robertson et al. Il2006allbl ). These gas-rich 
mergers also trigger star formation. Star forma- 
tion, in turn, adds dust to the interstellar medium; 
the dust attenuates star light. Therefore, in order 
to understand the Mbh-c^* relation for these ob- 
jects, as well as any other objects which are not 
dynamically relaxed, we need a better understand- 
ing of the evolution of cr* during the merger process 
as well as an understanding of how star formation 
and dust attenuation may influence this measure- 
ment. 

In the present paper, we have taken the first 
steps toward a theoretical understanding of the 
evolution of cr* during mergers of galaxies. We 



have also made initial steps toward understand- 
ing how dust attenuation might influence the mea- 
surement of fj* . Traditionally, simulation work in- 
volving Of has involved making meas urements of 
CTt af t er the merger is complet e (e.g., IGox et al~ 



2006t [Robertson et al. Il2006a llbh. We only know 



of one study in which ct* was plotted as a function 
of time during a merger ( Johansson et al.~ll2009l ) 
and, in that case, the measurement interval for 
cr* was 200 Myr, which (by the Nyquist sampling 
theorem) implies that the simulations could only 
capture the presence of fluctuations in ct* having a 
period of 400 Myr or more. To put this in context, 
we note that the dynamical timescales of the sys- 
tems in these simulations were less than 100 Myr. 
Furthermore, in simulations, cr* is often based on 
the velocities of the stars falling within the half- 
mass radius, rh, rather than on the velocities of 
stars appearing in a slit placed across the center 
of the system, as is the case in observational mea- 
su rements of a ^, (a n otable exception is the work 
of Cox et al. ( 2006[ ). which did employ a slit). 
The value of ct* that is typically reported is the 
mean value measured along 100-250 random lines 
of sight; the standard deviation of cr* over the set 
of viewing directions is typically not reported. Fi- 
nally, to our knowledge, no simulation study has 
attempted to measure the effect of dust attenua- 
tion on determinations of cr*. The work presented 
here is different in all of these regards. In our 
simulations, the mass- weighted and flux- weighted 
(T« were measured based upon the velocities of all 
stars appearing within a simulated diffraction slit 
centered on the nucleus of the galaxy. The flux- 
weighted measurement incorporated a toy model 
for dust attenuation. We measured cr* in this way 
along 10"^ lines of sight and performed a statisti- 
cal analysis on the directional distribution of cr*. 
We were particularly interested in identifying the 
time- variation of cr* during the merger process, 
thus these measurements were performed at very 
short time intervals during the merger. 

The paper is organized as follows. In Section 2, 
we describe the numerical simulations that were 
performed as well as the automated analysis al- 
gorithm that was used. In Section 3, we present 
the results of the numerical simulations and per- 
form some additional analysis. In Section 4, we 
conclude by discussing our findings, their implica- 
tions, and their limitations. 
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2. NUMERICAL SIMULATIONS 

In order to study the evolution of cr» with high 
time resolution during the merger process, we de- 
signed a simulation code that performs a statisti- 
cal analysis of the particle velocity data at short 
intervals during run time. We chose to keep the 
simulations as simple as possible so that we could 
identify the purely dynamical aspects of the evo- 
lution of cr,. The simulated galaxies were com- 
posed only of gravitationally bound, coUisionless 
star particles. There was no separate dark mat- 
ter component nor a gas component. At fixed 
intervals, the code computed mass-weighted and 
flux- weighted values of cr, along 10'^ random lines 
of sight. The measurements of tr* were based on 
the velocities of stars appearing in a rectangular 
"diffraction slit" centered on the projected center 
of mass of the system. This method allowed all 
stars along the line of sight to contribute to the 
measured value of oi — just as in the case of an 
observational measurement of ct*. Flux- weighted 
values of cr* were determined from the intrinsic 
luminosities of star particles, coupled with a toy 
model for dust attenuation. Using the values of 
(7* measured along 10"^ directions, the code com- 
puted the mean, minimum, maximum, standard 
deviation, skewness, and kurtosis of a,, during the 
merger process. 

2.1. Code Description 

Dynamics 

Our A^-body simulation code utilized the direct 
(i.e., all pairs) method to calculate the gravita- 
tional forces on each particle. The particles were 
softened in order to minimize the effects of two- 
body relaxation (White & Rees 1978). This en- 
sured that the stellar systems were coUisionless for 
the duration of the simulations. The acceleration 
of each particle, i, was computed using a truncated 
Plummer softening scheme, given by 



integrated forward in time using t he kick- drift 



(4 + 6^)-^/^ n,<e 



where r^j = — Vj and e is the softening length. 
Using this scheme, particles interact as Plum- 
mer spheres when overlapping significantly and 
as point particles otherwise. The system was 



kick form of the leap-frog integrator (jOuinn et al 



1997) with global, adaptive time stepping. 

Galaxy Construction 

In constructing progenitor galaxies, our code 
first initialized the stellar population. A list of 
stellar masses was made with relative abundances 
determined by an approxir nation to the Kroupa 
initial mass function (IMF: iKroupa 2001 ). Stars 
with masses M > 3Mq were then removed from 
the list in order to simulate the effect of an aged 
stellar population. Each of the N gravitationally 
softened particles in the simulation represented a 
population of n stars, thus the mass of each parti- 
cle was determined by randomly selecting n entries 
from the stellar population list. The luminosity 
of each particle was determined from the sum of 
the luminosities of the n constituent star particles, 
which were, in tur n, computed froni the mass- 



luminosity relation ( Salaris fc Cassisi 20061) 



Once the particles were initialized, they were 
distributed in space according to the Hernquist 



density profile ( Hernquist 1990), 



p{r) 



M 



D 



2Tir{r + Df 



(1) 



where M is the total mass of the system, and the 
size scale, £), of the distribution is the "Hernquist 
radius." The velocity distribution was chosen to 
reproduce the analytic result obtained from the 
Jeans equations, assuming an irrotational system 
with an isotropic dispersion tensor and Hernquist 
density profile (see the Appendix). The resulting 
system, having been constructed stochastically, 
exhibited a small net linear momentum, implying 
a small drift velocity. All particle velocities were 
adjusted to remove the net linear momentum. Our 
newly constructed system was then evolved for- 
ward in time. The density profile and central ve- 
locity dispersion of the newly constructed system 
varied for approximately 2.5 dynamical timescales 
before reaching a steady configuration. This pe- 
riod of adjustment was apparently necessitated by 
the mismatch between the smooth analytic distri- 
bution used to initialize the system and the gran- 
ular approximation that was actually produced. 
The final density profile of the progenitor system 
fit the Hernquist profile perfectly — within the lim- 
its of particle noise. After the initial 2.5 dynamical 
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timescales, no evolution was observed during the 
subsequent « 290 dynamical times for which the 
system was studied. Based on this analysis, we 
evolved each progenitor spheroid passively for 3.0 
dynamical times before using it in a merger simu- 
lation. 

Dust Attenuation 

In order to determine how measurements of tr* 
might depend on dust attenuation, a toy model 
for dust attenuation was included in the code. A 
cylindrical slab of gray attenuating material of ra- 
dius, r, thickness, (5, and extinction coefficient, n 
was placed concentric with the center of mass of 
the merging system. The flux of the «th star in 
the system, /j, was then calculated according to 



fi = Lie 



— hidi 



(2) 



where Li is the luminosity of the ith star particle 
and di is the distance that light from the ith star 
traveled through the attenuating slab on its way 
to the virtual observer. 

Measuring cr* 

The process used to measure ct* in our code re- 
sembled the common observational technique, in 
that a rectangular slit of width w and length £ 
was centered on the system and the measurement 
of (7* was based on the line-of-sight components 
of the velocities of stars within the slit. This 
method is different from the one often employed 
in many numerical simulations; simulators typi- 
cally base measurements on the velocities of 
stars within the halfmass radius or projected half- 
mass radius of the system ( Johansson et al. Il2009t 



Robertson et al. Il2006al lbl). 

The measurement of cr* proceeded as follows. 
First, a random viewing direction (0, </>) was cho- 
sen, then the system was rotated such that the new 
+z-axis coincided with the {0,(j)) direction. The 
set of particles within the observing slit was identi- 
fied using the new x- and ^-coordinates. The line- 
of-sight velocity of each particle was then given 
by the new z-component of velocity. From the 
line-of-sight velocity distribution within the slit, 
the mass-weighted and flux-weighted versions of 
cr* were computed as follows 

mcr* = ^Jvfm.,/M ~ {v^m^/Mf (3) 



with 



(4) 



where the standard summation convention has 
been utilized; repeated indices imply a sum over 
that index. We use the notation mtTn, and fa* to de- 
note mass- weighted and flux- weighted cr*, respec- 
tively. This process was repeated for 10^ random 
directions. 

Directional Statistics 

Once mcr* and fcr* were measured for 10^ di- 
rections, statistical quantities were computed in 
order to determine the degree of anisotropy of the 
merger system. Speciflcally, the mean, minimum, 
maximum, standard deviation, skewness, and kur- 
tosis of mcr* and fcr* were computed for the set of 
directions. Since the deflnitions of skewness and 
kurtosis vary among authors, we present the defl- 
nitions that we used below 

skew(x) = n"^ ^(xi - ^)^/cr^ (5) 

i 

kurt(x) = n~^Y.^x,- txf/a^ (6) 

i 

where ^ and a are the mean and standard devia- 
tion of a;, respectively and n is the sample size. 
The directional statistics, combined with time- 
evolution data, enabled us to estimate the intrinsic 
scatter in observations of a* for randomly oriented 
active mergers and merger remnants. 

Collisions 

The code was specifically designed to study the 
simplest type of mergers: binary mergers of equal- 
mass systems (i.e., 1:1 mergers). Once a model 
galaxy was constructed, as described above, the 
galaxy was replicated and placed on a collision 
course with its clone. The relative position and 
velocity vectors (R and V) were specified, then 
the code adopted the zero-linear-momentum ref- 
erence frame and set the origin of the coordinate 
system coincident with the center of mass. This 
ensured that the final merger remnant would be 
centered at the origin. 
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2.2. Precision 

Particle noise was by far tlie main source of 
uncertainty in our measurement of cr* . We quanti- 
fied the precision in the primary measurement of 
interest — the directional mean of cr* — through ex- 
perimentation. First, we noted that our newly 
constructed galaxies were perfectly spherically 
symmetric and isotropic — except for the statis- 
tical noise introduced by using a finite number 
of particles, N. In the limit as iV — > cxi, the 
measurement of cr* should be independent of the 
direction from which the measurement was made. 
This implies that the standard deviation of cr* 
over the set of all viewing directions (tJd) should 
vanish in this limit. Based on this insight, we 
performed measurements of the standard devia- 
tion of cr* in newly constructed systems, measured 
from 10"^ random directions, for various N ranging 
from 10^ to 10^. We found the expected behav- 
ior: (Td N^^/"^. Upon determining the constant 
of proportionality associated with our simulation 
parameters (namely, the system density and slit 
dimensions), we were able to quantify the degree 
of particle noise in our measurements of cr, . 

The uncertainty introduced by the error in 
the numerical integration scheme was negligible 
in comparison with the particle noise. All non- 
integer quantities in the code were stored as dou- 
ble precision floating point numbers. We used a 
very stringent step-size criterion in which no par- 
ticle was allowed to move more than l.le during 
a time-step. The total momentum of the system 
was conserved to within a small multiple of the 
machine precision and the total energy fluctuated 
by less than 0.1%. 

2.3. Simulation Units & Parameters 

Internally, the code used a system of units in 
which mass was measured in solar masses, dis- 
tance was measured in parsecs, the gravitational 
constant G = 1, and time was a derived unit. 
However, in order to make our results easier to in- 
terpret, we have re-scaled the simulation param- 
eters. All quantities will be presented in terms 
of the characteristics of the progenitor galaxies. 
Our unit of distance is the half- mass radius, rh, 
of the progenitor, which can be written in terms 
of the Hernquist radius, D, as rh = (1 + ^/2)D. 
For a system with constant mass-to- light ratio, the 



projected ha lf-light radius, R p, is related to rh by 
Th « 1.33i?o (|Hernquist ][l993). Our unit of time is 
the dynamical time-scale of the progenitor, given 

by 



1 ^dyn 



'7^2(l + ^/2)3j:)3 
2GM 



(7) 



16Gph V 

where ph is the mean density within rh. The unit 
of velocity is then given by 



1 



2GM 



7r2(l -f \/2)D 



(8) 



where M is the total mass of the system. In terms 
of these units, the gravitational constant is. 



G 



2 Mt%^ 



(9) 



To gain a better understanding of this unit system, 
consider a spheroid of mass 1.0 x 10^" M© with an 
effective radius of i?c — 2.0 kpc. For this spheroid, 
the time unit is 37 Myr and the velocity unit is 
70 km s~^. 

After studying a large variety of initial orbital 
parameters, we found that the general behavior of 
all of the merger simulations fell between two ex- 
treme cases: head-on collisions and orbital decay 
mergers. In a head-on collision, the progenitors 
are initially separated by some distance, i?, and 
given an initial relative speed, V , directly toward 
one another. In the orbital decay scenario, the two 
galaxies begin on circular or nearly circular orbits. 
As they interact through tidal forces and dynam- 
ical friction, angular momentum is redistributed 
and the galaxies gently merge. We analyzed these 
two extreme cases, as well as a representative in- 
termediate case, in detail. For all simulations, 
the radius, thickness, and attenuation coefficient 
of the cylindrical attenuating slab were r = 3.36 
Th, 6 = 1.24 Th, and n = 1.8 Tj^^, respectively. 
The slit width and length were w = 0.14 rh and 
£ = 0.70 rh, respectively. Table [1] summarizes the 
parameters that were varied in our simulations. 

3. RESULTS 

3.1. Measurement Techniques 

Slit (T* versus Half-mass cr* 

Using our slit-based method of measuring cr*, 
with the slit dimensions given in Section 12.31 
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Table 1: Summary of Parameters Used 



Simulation Type 


N 






R (rh) 


9 (deg) 


T (tdyn) 


Head-on, short 


50,000 


0.0233 


1.55 


1.66 


180 


0.00942 


Head-on, long 


35,000 


0.0262 


1.55 


1.66 


180 


0.0942 


Intermediate, short 


50,000 


0.0233 


1.60 


1.69 


153.2 


0.00942 


Intermediate, long 


35,000 


0.0262 


1.60 


1.69 


153.2 


0.0942 


Orbital decay, short 


50,000 


0.0233 


1.55 


1.66 


90 


0.00942 


Orbital decay, long 


35,000 


0.0262 


1.55 


1.66 


90 


0.0942 



Note. — The qualifiers "short" and "long" refer to the duration of the simulation and the interval between measurements of 
cr*. All "short" simulations were evolved forward for fSJ 8 tdyn while the "long" were allowed to evolve for 290 tdyn- The 
number TV is the number of particles in each progenitor system; the total number of particles in the merger simulation is 2N. 
The initial relative speed of the galaxies is V = |V|, their initial separation distance is R = |R|, and the angle between R and 
V is 0. The interval between measurements of cr, is t. 



we found that cr* in our progenitor systems was 
cr* — 1.07 ± 0.01 Tht^yjj while the value of cr* 
measured within the central spherical region of 
radius rh was 1.0047 ± 0.0026 rh t'^^. For com- 
parison, the result obtained by analytically com- 
puting the mass- weighted mean of the line-of-sight 
velocity dispersion within a sphere of radius rh for 
an isotropic Hernquist profile is ~ 1.0035 rhijy„ 
while the same quantity computed within rh/2 is 
~ 1.0693 rh ijyjj (see the Appendix). Evidently, 
measuring a* using a narrow slit placed on the 
center of the galaxy yielded values that were closer 
to the velocity dispersion within a region smaller 
than rh. As the width of the slit was increased, 
the difference between the two measurement tech- 
niques diminished. In quiescent, isotropic systems, 
the two measurement techniques agreed within 
2% when the width of the slit it; = rh and the 
methods agreed to within the measurement un- 
certainty when w = 2rh. As the width of the 
slit approached zero, particle noise made the mea- 
surement highly unreliable. For non-quiescent and 
non-isotropic systems, the two methods sometimes 
differed significantly — even when the slit was quite 
large. This disagreement was due to the motion 
of particles lying outside of the spherical half-mass 
region, which were included in the slit-based mea- 
surement, but absent from the spherical half- mass 
measurement. 

Flux Weighting vs. Mass Weighting 

The presence of an attenuating slab placed con- 
centric with the galaxy had the effect of reducing 



the flux received from stars in the central region 
of the galaxy. Using the properties of the atten- 
uating slab presented in Section 12.31 along with 
Equation we see that the flux received from 
a star particle at the center of the galaxy was 
/ = 0.3275L when viewed through the thinnest 
portion of the slab (face-on), where L is the in- 
trinsic luminosity of the star particle. When 
viewed through the thickest portion of the slab 
(edge-on), / = 0.0025L. When velocity disper- 
sions were computed using the flux- weighting tech- 
nique, the velocities of stars in the central region 
of the galaxy were therefore weighted less heav- 
ily. Since the stellar velocity dispersion was largest 
near the central region of our simulated galaxies, 
the flux- weighted measurement, fcr*, was smaller 
than its mass- weighted counterpart, mcr*. Using 
a simulation with 10^ particles, we found that 
fcr* = (0.875 ± 0.015)m(T*. In general, increasing 
the attenuation coefficient or the dimensions of the 
attenuating slab caused fcr* to decrease monoton- 
ically. 

The presence of the cylindrical attenuator also 
destroyed the spherical symmetry of the system; 
when the mass-weighted velocity distribution was 
isotropic, the flux-weighted velocity distribution 
was not. As mentioned in Section [521 the stan- 
dard deviation of cr* over the set of random view- 
ing directions (cTd) can be used as a measure of 
the anisotropy of the system's velocity distribu- 
tion. Using the same 10^ particle simulation men- 
tioned above, the standard deviation of fcr* over 
the set of 10'^ random viewing directions was fcTd = 
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(0.123 ± 0.017)(T* while rndd = (0.006 ± 0.001)cr*. 

Additionally, we found that changing the 
radius-to-thickness ratio ol the attenuating slab 
{r/5) caused the anisotropy of the flux- weighted 
velocity distribution, fcd to change. The minimum 
value of ffJd was found when r/5 ~ 0.5, of course 
f(Td also tended toward zero if either — > or 
r/rh — >■ 0, regardless of the ratio r/5. The value 
ol r/5 leading to maximum isotropy occurred in 
the interval \ < r / 5 < oo] the exact ratio r/5 that 
maximized fcrd depended on the relative size of 
the slab compared with that of the galaxy. 

3.2. Merger Evolution 

The dissipationless mergers that we studied 
proceeded in three primary stages: oscillation 
then phase mixing, and finally completion or dy- 
namical equilibrium. Figures [TH3] respectively il- 
lustrate the early stages of the "Head-on, short" , 
"intermediate, short", and "orbit-decay, short" 
simulations from Table [1] 

Oscillation 

The oscillatory stage is characterized by the 
bulk motion of the two progenitor galaxies as they 
coalesce. As the progenitor nuclei become super- 
imposed significantly for the first time, the density, 
gravitational potential, and tr* increase rapidly. 
The progenitor nuclei then typically pass through 
one another. The density, gravitational potential, 
and oi of the central region of the merging system 
are then temporarily reduced. The nuclei of the 
systems eventually change directions and fall back 
onto one another, while some of the star parti- 
cles that were initially less tightly bound continue 
on their original paths — only slightly perturbed by 
the motion of the nuclei. The process then repeats 
several times until the coherent oscillations decay 
away. The evolution of cr* then becomes domi- 
nated by phase mixing. In the Head-on merger 
shown in Figure [1] and the Intermediate merger 
shown in Figure [21 the transition between the os- 
cillatory and phase mixing stages occurred at ap- 
proximately 1.5 tdyn and 2.3 tdyn, respectively. In 
the less violent Orbital decay merger shown in 
Figure [3l the oscillatory stage was much less pro- 
nounced because the galaxies gently spiraled into 
one another, but oscillations were still visible be- 
tween 2.0 tdyn and 4.0 tdyn- 



Phase Mixing 

The phenomenon of phase mixing arises be- 
cause of the dispersion in the oscillation (or or- 
bital) periods of the particles composing a system. 
Due to this dispersion, an initially coherent oscilla- 
tion becomes incoherent with time and eventually 
vanishes. For a system consisting of a finite num- 
ber of particles, this sort of phase mixing alone 
cannot be fully responsible for the termination of 
the oscillation stage described above because the 
differences between the orbital periods are neces- 
sarily finite. A finite system of uncoupled, un- 
damped, undriven oscillators will periodically re- 
turn to its initial state, meaning that episodes 
of coherent oscillation would recur periodically. 
Phase mixing plays an important role in ending 
the oscillation stage, but particle interactions are 
essential in preventing the recurrence of coherent 
oscillations. The stars in a galaxy interact in- 
directly through a time-varying global potential. 
This sort of interaction, which is most obvious 
during the oscillatory s tage of a merger, is known 
as "violent relaxation" (|Lvnden-Bell 1119671 ) . Stars 
also interact through close encounters with other 
stars, however this mode of interaction is insignif- 
icant compared with violent relaxation for most 
stars in a typical galaxy during the timescale of 
a merger. These interactions cause the period of 
each particle to vary with time in a non-periodic 
way which modifies the pure phase-mixing mech- 
anism and prevents the return of coherent oscilla- 
tions. 

Visually, the process of phase mixing causes the 
particles of the two progenitor systems in a merger 
to lose their identities; it gradually becomes im- 
possible to distinguish the distribution of particles 
that initially belonged to progenitor A from the 
distribution of particles that originally belonged 
to progenitor B. T his is illustrated pa rticularly 
clearly in figure 5 of Funato et al. I (|l99 2a'). 

In our Head-on and Intermediate merger simu- 
lations we observed many statistically significant 
variations in ct* during the phase mixing stage 
of evolution. These fiuctuations can be seen in 
Figures [U [51 S] and [SI We define the end of 
the phase mixing stage as the time after which 
the directional mean of cr* can be found within 
the la particle noise of the equilibrium value, 
c*, final ± Cd, noise, with 75% Confidence (note that 
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time (idyn) 

Fig. 1. — Evolution of cr, with time during the "Head-on, short" merger of Table [1] We have scaled a, by the 
value of velocity dispersion of the progenitor galaxies (cTprog)- The black curve represents the directional mean of the 
mass- weighted line-of-sight velocity dispersion (mo-*) measured using the slit-based method. The upper and lower 
limits of the gray region are, respectively, the maximum and minimum values of ma, for the set of 10^ random 
viewing directions. The gray horizontal stripe is centered on the final equilibrium value of mcr, of the merger 
remnant. The half-thickness of the stripe illustrates the la uncertainty due to particle noise. Thus the gray stripe is 
(ma* final ± crd,noisc)/fpi.og ~ 1.440 ± 0.008. The centers of the two progenitor galaxies first coincide at t « 0.5 tdyn. 
The oscillation stage ends at t ~ 1.5 tdyn- The velocity dispersion continues to fluctuate, although less significantly, 
in the phase mixing stage. Note: Although this simulation was evolved forward for ~ 8 tdyn, the plot only shows the 
flrst 4 tdyn in order to highlight the oscillatory stage of the merger. 
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time (idyn) 

Fig. 2. — Evolution of ct* witii time during the "Intermediate, short" merger of Table [1] The black line, shaded 
region, and gray horizontal stripe are defined the same as in Figure [T] The position and thickness of the gray stripe 
are (ma.^flnai ± o"d,noisc)/o"prog = 1.373 ± 0.008. The centers of the two progenitor galaxies first coincide at t ~ 0.7 
tdyn- The oscillation stage ends at t « 2.3 tdyn- 




0.8 - 



0.6 



time (idyn) 



Fig. 3. — Evolution of a, with time during the "Orbit decay, short" merger of Table[T] The black line, shaded region, 
and gray horizontal stripe are defined the same as in Figure [T] The position and thickness of the gray stripe are 
(mcr,_finai ± o"d,noiso)/o"prog = 1.361 ± 0.008. The centers of the two progenitor galaxies first coincide at t ~ 3.1 tdyn- 
The oscillation stage ends at t « 4.0 tdyn- 
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the process of phase mixing continues to occur 
ad infinitum — only the phase mixing stage has an 
end) . This is admittedly a somewhat arbitrary cri- 
terion. For the Head-on collision shown in Figures 
m the phase mixing stage ended at « 17tdyn- In 
the Intermediate simulation, shown in Figures [Sj 
the end of the stage occurred at « llidyn- The 
Orbital decay simulation, shown in Figure [3] did 
not have a clearly distinct phase mixing stage. In 
Figures m and [S] a small number of 2(t fluctuations 
in cr* can be seen after the stated end of the phase 
mixing stage. Relatively large fluctuations, such 
as these, occurred very infrequently after the end 
of the stage compared with their frequency before 
the end of the stage. 

During the phase mixing stage, fluctuations in 
the mean value of ct* were accompanied by corre- 
sponding fluctuations in the minimum and maxi- 
mum values of a^,. This finding is consist ent with 
the work of Merrall k. Henriksen ] (l2003l) . which 
found that the virial ratios, 2T/W oscillated with 
time as the systems progressed toward dynamical 
equilibrium (T and W are the kinetic and poten- 
tial energy, respe ctively) . Simila r oscilla tions were 
also discussed bv lFunato et al. ( 1992al b). These 
oscillations appear to be due to the presence of 
small subsets of particles in the merging system 
with very similar orbital periods. The rate of 
phase mixing depends monotonically on the dif- 
ference between periods, thus particles with very 
similar periods mix slowly, which explains why 
these smaller oscillations survive longer than the 
bulk oscillation s in th e in itial oscillatory stage 



(iFunato et al. I Il992bh . iMerrall fc Henriksen 



2003 ^ found that such oscillations diminished in 
the continuum limit, as — )■ oo, but they noted 
that the fluctuations are indeed statistically sig- 
nificant and represent a physical phenomenon that 
occurs in systems consisting of a finite number of 
particles. 

Note that the only differences between our 
"short" and "long" simulations (from Table [T]) 
were the number of particles, the gravitational 
softening length, and the measurement inter- 
val, r. The initial orbital parameters, total mass, 
size, and shape of the density profile were iden- 
tical. This enabled us to determine whether the 
fluctuations during the phase mixing stage differed 
for systems having the same macrostate, but dif- 
ferent microstates. We compared the evolution of 



(J* during the interval shared by both the "long" 
and "short" simulations (i.e., the flrst 8 idyn of 
each simulation) and found that the variations in 
cr* were in perfect agreement within the noise lim- 
its and sampling frequency. Notably, every fluc- 
tuation in the "long" simulation was accompanied 
by a corresponding fluctuation in the "short" sim- 
ulation. Of course not every fluctuation present 
in the "short" simulations could be detected in 
the "long" simulations because the sampling fre- 
quency in the "long" simulations was a factor of 10 
lower and the noise threshold was slightly higher. 
Furthermore, the pseudorandom number gener- 
ator function that was used in the galaxy con- 
struction algorithm was initialized with a differ- 
ent random seed each time a new spheroid was 
constructed. This caused the particles to be ar- 
ranged differently each time a new system was 
constructed. Therefore, the similarity between the 
"short" and "long" simulations implies that the 
statistically signiflcant fluctuations observed dur- 
ing the first 8 tdyn of the "long" simulations were 
evidently not artifacts of the detailed microstates 
of the particle systems. The fluctuations depended 
on the macrostate — the density profile, velocity 
distribution, and initial orbital parameters. 

As noted previously, the phase mixing stage in 
the Head-on collision lasted longer than the phase 
mixing stage in the Intermediate collision. The 
Orbital decay merger lacked a distinct phase mix- 
ing stage. It appears that the duration of the 
phase mixing stage is inversely related to the dura- 
tion of the oscillation stage. The magnitude of the 
time derivative of the mean gravitational potential 
in the early stages of the merger may be responsi- 
ble for this relationship, however much more anal- 
ysis would be required to verify this hypothesis. 
The efficiency of the mixing process may also play 
a role in determining the duration of the phase 
mixing stage. It is clear that the merger with 
the largest angular momentum (Orbital decay) ex- 
hibited the most short-lived phase-mixing stage, 
while the merger with zero angular momentum 
(Head-on) exhibited the longest-lived phase mix- 
ing stage. The (differential) rotation due to the 
presence of angular momentum, combined with 
dynamical friction may have accelerated the mix- 
ing process. 
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Fig. 4. — The evolution of a, witli time during the "Head-on, long " merger of Table[T] The black line, shaded region, 
and gray horizontal stripe are defined the same as in Figure [T] The position and thickness of the gray stripe are 
(ma",,flnai ± o-d.noisc)/o"prog = 1.4399 ± 0.0096. The phase mixing stage ends at f 17 tdyn- The vertical dotted lines 
mark the approximate boundaries between the oscillation, phase mixing, and dynamic equilibrium stages. Note that 
the maximum value of a* (the upper boundary of the shaded region) continues to evolve until « 45 tdyn. Although 
this simulation was evolved forward for ~ 290 tdyn, the plot only shows the first 60 fdyn because there was no visible 
evolution beyond 60 tdyn- 




time (idyn) 



Fig. 5. — The evolution of cr* with time during the "Intermediate, long " merger of Table [T] The black line, shaded 
region, and gray horizontal stripe are defined the same as in Figure \T\ The position and thickness of the gray stripe 
are (mcr«,finai ± o"d,noisc)/crprog ~ 1.3730 ± 0.0096. The phase mixing stage ends at t 11 tdyn. 
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Dynamical Equilibrium 

Once the phase mixing stage was complete, 
we observed no significant evohition in the mean 
value of cr* in the remnant systems during the 
course of the « 290 tdyn that the simulations 
were allowed to evolve. We determined the re- 
laxed values of the mass-weighted velocity disper- 
sion (mcr* final) by computing the time average 
of tr* during the final « 200 tdyn of the "long" 
simulations. The mean value remained approx- 
imately constant — only rarely exceeding the la 
noise limit. For the Head-on merger, mcr*_finai = 
1.54 ± 0.01 rhtjy„. The Intermediate and Or- 
bital decay mergers had identical velocity disper- 
sions within the measurement uncertainty, with 
m(T*,fi„ai = 1.469 ± 0.010 Thtj/n and mcr^^nai = 
1.456 ± 0.010 rh tjy„ respectively. As expected 
from the initial conditions, the Head-on merger 
remnant exhibited no net rotation; it was entirely 
supported by pressure. The Intermediate and Or- 
bital decay remnants were supported by both ro- 
tation and pressure. 

We observed t he same general remnant proper 



ties de scrib ed bvlGonzalez-Garcia fc van Albada 



(|2005[ ) and IVillumsen I (|1982| ). In particular, for 
the Head-on collision remnant, the velocity dis- 
persion was largest along the collision axis and 
smallest perpendicular to the axis. For the ro- 
tating remnants, the mean velocity dispersion in 
the orbital plane of the collision was larger than 
the velocity dispersion perpendicular to the orbital 
plane. 

Figures [BHH] show the evolution of mcr* and fu* 
along with their corresponding directional stan- 
dard deviations mad and fcTd for our three merger 
simulations. While there was no evolution in the 
mean of ct* , the directional distribution of oi and 
the maximum value of cr* over the set of directions 
did evolve somewhat during the dynamical equilib- 
rium stage. This evolution is particularly evident 
in Figures HI [6l and [7l The anisotropy (ud) and 
the maximum value of cr, decreased and eventu- 
ally reached stable values a few tens of dynamical 
timescales after the end of the phase mixing stage. 

3.3. Additional Statistics 

Suppose a non-quiescent galaxy is observed and 
a careful measurement of cr* is made. With what 
probability will the measured cr* fall within a spec- 



ified range ±A(t» of the system's eventual equilib- 
rium value of cr*, (i.e., the value of cr* nominally 
used in the Mbh-ct* and FP relation studies)? In 
other words, what degree of scatter is expected 
when measuring cr* in a non-quiescent system and 
is there an offset in cr* during dynamically non- 
quiescent times? The results presented above pro- 
vide the first pieces of information needed to an- 
swer these questions. For instance, in Figure [I] 
we see that, even in the most violent of the 1:1 
mergers that we studied, the maximum value of 
cr, that would ever be measured is less than a fac- 
tor of two greater than the final equilibrium value. 
Neglecting the effect of flux- weighting, the mini- 
mum value of cr* measured in an apparently coa- 
lesced system was about 70% of the equilibrium 
value. We expect that cr* varies most significantly 
in 1:1 mergers, therefore this work suggests that a 
measurement of cr* made in a non-quiescent galaxy 
would fall between 70% and 200% of the quiescent 
value. Furthermore, a random measurement of cr* 
is much more likely to fall near the equilibrium 
value than near an extremum; in order to obtain 
a measurement of tr* near an extremum, the sys- 
tem would have to be observed from a fortuitous 
viewing angle during a fairly short epoch. 

The standard deviation of a^, for the set of view- 
ing angles in each merger is shown in the lower 
panels of Figures EHU In these plots, we see that — 
in all three simulations — mad was largest during 
the oscillatory stage with an absolute maximum 
value of about 25% of the equilibrium value of 
cr*. After the oscillatory stage was complete, ma^ 
was typically less than 3% of the equilibrium value 
of Of. The flux- weighted quantity, fcd was larger 
than mcd due to the anisotropy introduced by the 
slab of attenuating material. 

We have not presented plots of the skewness or 
kurtosis of the directional distribution of fT^ as a 
function of time because the large amount of noise 
present in these quantities prevented us from de- 
tecting any trends with time. There was, how- 
ever, a clear difference between the mass- weighted 
and flux- weighted measurements. These differ- 
ences are summarized below. 

The skewness of the directional distribution of 
mcr,, in the isolated progenitor system was 0.01 ± 
0.23, which is consistent with zero, as expected 
from the manifest symmetry of the progenitor. Its 
flux-weighted counterpart was —2.8 ± 0.2. The 
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time (tdyn) 



Fig. 6. — The evolution of a, with time during the "Head-on, long " merger of Table [T] Upper panel: The black 
and gray lines show the evolution of moi and fcr* with time, respectively. The flux-weighted quantity is systematically 
lower than its mass-weighted counterpart. Lower panel: The black and gray lines show the dispersion in the mass- 
weighted and flux- weighted values of a, for a set of 10'^ random viewing directions (crd). This is effectively a measure 
of the anisotropy of the system. The lower dotted line shows the particle noise threshold which represents a lower limit 
on (Td5 a perfectly isotropic system composed of 7 x 10^ particles would have a value of f7"d/ci"prog = 0.0096, which is the 
position of this line. From this plot, we can see that flux-based measurements would indicate that the system is much 
less isotropic than it actually is. This means that flux-based measurements of ct, from different viewing directions 
could differ significantly. The mass-weighted measurement shows that the system continued evolving toward a more 
isotropic state well after the end of the phase mixing stage. 
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Fig. 7. — The evolution of ct, witii time during tlie "Intermediate, long " merger of Table[T] See the caption of Figure 
|6]for a description of the plotted quantities. The effects of flux- weighting are the same as in the Head-on merger. 




time (tdyn) 



Fig. 8. — The evolution of cr* with time during the "Orbital decay, long " merger of Table [T] See the caption of 
Figure [6] for a description of the plotted quantities. Note that the remnant system shown here is less isotropic than 
the Head-on and Intermediate merger remnants. This is because the system is somewhat flattened due to rotation. 
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skewness of mcr* was 0.49 ± 0.23 for each of the 
three simulations during the oscihatory and early 
phase mixing stages of evolution, but the final qui- 
escent values differed somewhat. The quiescent 
values of skewness for the Head-on, Intermediate, 
and Orbital decay merger remnants were respec- 
tively 0.07±0.22, -0.37±0.25, and -0.53±0.23. In 
each simulation, the skewness of the flux- weighted 
distribution was — 1.0 ± 0.2 during the oscillatory 
and early phase mixing stages and — 3.1 ± 0.1 for 
each of the remnant systems — consistent with the 
skewness of fcr* in the progenitors. Thus, the 
skewness of the flux-weighted velocity dispersion 
depended more strongly upon the dust geometry 
than the actual stellar dynamics. 

Merger evolution had no measurable effect on 
the kurtosis of the directional distribution of mcr* . 
The kurtosis of mcr* was 2.5 ± 0.4 for the progen- 
itors and remained fixed at this value, within the 
limits of noise, for the duration of all three merger 
simulations. The kurtosis of fa^ in the progenitors, 
as well in each merger simulation, was consistently 
11.6 ± 1.9. This elevated kurtosis indicates that 
much of the variability in the directional distribu- 
tion of foj was due to extreme outliers, which is 
a direct result of the geometry of the attenuating 
slab. 

4. DISCUSSION AND CONCLUSIONS 

By analyzing the evolution of cr* during three 
coUisionless, dissipationless, 1:1 mergers of spher- 
ically symmetric systems, we identified three pri- 
mary stages of evolution. During the most dra- 
matic, early stage of the merger, cr* undergoes 
large damped oscillations of increasing frequency. 
Following the oscillation stage, the value of cr* 
can fluctuate significantly in an apparently chaotic 
way for more than 10 dynamical timescales as the 
system becomes more mixed. We called this the 
phase mixing stage. The phase mixing stage is 
followed by the dynamical equilibrium stage dur- 
ing which the value of cr* remains essentially fixed, 
while the system evolves toward a final equilib- 
rium state. Using the statistics computed during 
the collisions, we identified the extreme limits of 
cr* during such galaxy mergers and provided esti- 
mates of the scatter inherent in making measure- 
ments of fj* at random times from rand om viewing 
angles during a merger. The work of ICox et al 



( 20061 ) hints that the evolution of tr* may be more 
complicated in mergers that contain a dissipative 
component (i.e., gas). If a similar analysis were 
performed on a large variety of more realistic dis- 
sipative galaxy mergers with varying mass ratios, 
gas fractions, Hubble types, and initial orbital pa- 
rameters, it would be possible to predict the scat- 
ter inherent in making measurements of cr* in non- 
quiescent systems in general. 

By measuring cr* using a slit-based measure- 
ment method coupled with a toy model of dust 
attenuation, we found that the presence of dust 
in a galaxy can systematically decrease the fiux- 
weighted value of cr* relative to the mass- weighted 
measurement. Furthermore, the distribution of 
dust can increase the apparent anisotropy of a sys- 
tem. This increases the observed scatter in obser- 
vational determinations of (T, . In order to under- 
stand how dust infiuences measurements of in 
real galaxies, a more realistic model for dust at- 
tenuation is needed. We have begun a follow-up 
project th at will use the radiative transfer c ode 
Sunrise, (| Jonsson II2006I : Ijonsson et al. Il20l d) to 



create Doppler- broadened spectra of GADGET-2 
(jSpringel 1120051) merger simulations at fixed time 
intervals. The synthetic spectra will be analyzed 
to determine the fiux- weighted value of tr*. This 
will allow us to compute ct* in a way that is fully 
consistent with the method used by observers; not 
only will the cr* measurement be fiux-based, but 
cr* will be obtained by fitting spectral line profiles 
as opposed to using direct particle data as in the 
present work. In addition to providing us with 
a better understanding of how dust effects mea- 
surements of cr*, the follow-up project should also 
allow us investigate how star formation infiuences 
the measurement of tr*. Rothberg & Fischer re- 
cently reported a systematic discrepancy between 
values of cr* measured at different wavelengths, 
with lower measured fo r longer wavelengths 
(|Rothberg fc Fischer I [ioioh . Their proposed ex- 
planation for this "cr* discrepancy" is as follows: 
stars for which cr* was measured using near-IR 
CO lines are young stars located in a dusty ro- 
tating gaseous disk whereas stars measured using 
the shorter wavelength Ca ill triplet were older. 
The young stars, having recently formed from col- 
lisional, dissipational, molecular gas, have lower 
velocity dispersion than the surrounding popula- 
tion of old stars because the clouds from which 
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they formed had lower velocity dispersion than 
the older stellar population. These stars have not 
yet had time to mix with the older population 
and adopt the higher cr*. This expl anation seems 



number AST 0507450. 



consistent with the observations of Genzel et al. 



()200lh which found that the gas dynamics and 
stellar dynamics become decoupled during merg- 
ers. Our follow-up work should aid in understand- 
ing this "sigma discrepancy" . 

The merger simulations described in this work 
were intentionally kept simple in order to allow 
us to identify the most fundamental, purely dy- 
namical aspects of the evolution of during a 
merger. The systems were spherical, isotropic, 
non-rotating, and contained no gas nor dark mat- 
ter. In spite of the simplicity, we observed nontriv- 
ial aspects of the evolution. In order to identify 
how each additional bit of complexity effects the 
evolution of cr* , we will compare these fundamen- 
tal aspects of the evolution with the results of the 
more realistic simulations in our follow-up project. 
Even without performing further simulations, we 
can see that the time-scales and the set of pos- 
sible initial orbital parameters will increase when 
we add a dark matter component. This is because 
systems of stars embedded in large dark matter 
halos are able to eventually merge as long as their 
parent halos interact sufficiently; the stellar sys- 
tems themselves do not need to become super- 
imposed during the earliest stages of the merger 
process because the earliest interactions primarily 
take place in the outer regions of the dark matter 
halos. The systems we modeled needed to start out 
on a collision course in order to merge. Therefore, 
our simulations only represented the later stages 
of the overall merger process — beginning with the 
stage at which a collision between the stellar com- 
ponents was imminent. 
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A. Analytic Results 

The mass-weighted mean value of the radial velocity dispersion in a region of volume V is given by 

mcr* = -v/mof (Al) 

with 

where My is the mass within the volume V, p is the mass density, and (u^) = ((f • Vcm)^) is the mean 
squared radial component of the velocity with respect to the center of mass of the system (i.e., the velocity 
in the zero momentum frame) . For a spherically symmetric system, this can be written in terms of r as 



mcr^(r) 



47r 
M(r) 



{v'',)p{r')r'^dr' 



(A3) 



where M(r) is the mass enclosed within radius r. A system having an isotropic velocity dispersion tensor 
throughout will, by definition, have the same velocity dispersion in all directions at any arbitrary position 
r within the system. Thus, for an isotropic system, we could measure the velocity dispersion along any 
direction, compute the mass- weighted mean, and arrive at the same result for mcr*. In particular, if we 



choose the line-of-sight direction, n, then (w„) = ((n • Vrm) ) = ((i 



an isotropic, non-rotating Hernquist profile is given in iHernquist (1990) and reproduced below 



) = (Wj ). An expression for {v^) in 



GM f 12r(r + Df ^^^(r + D 



12D 



D 



r / r 

25 + 52— + 42 ( — 



D 



Adopting a system of units in which G = D = M = 1, 



In these units, 1 rh = 1 



(vf) = r{r + if 
\/2 and 



ln(l 



^dyn 



4]) becomes 
25 + 52r + 42r2 



\D 



12H 



12(r + 1)4 



7r2(l + y2) 



(A4) 



(A5) 



(A6) 



Upon substituting (jA5P into (|A3|) with r = rh, evaluating the integral, computing the square root, and using 
l6p to express the result in terms of rh and idyn, we arrive at the result. 



mcr* w 1.0035 



idyn 



Substituting r = rh/2 yields 



mcr, w 1.0693 



idyn 



(A7) 
(A8) 
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